function [Gamma, F, iter] = IPCA(X,W,Nt,K,varargin)

% Standard IPCA call as IPCA(X,W,Nt,K)
% IPCA with alpha call as IPCA(X,W,Nt,K,ones(1,T))
% IPCA with PSF call as IPCA(X,W,Nt,K,PSF)
% IPCA with ONLY PSF call as IPCA(X,W,Nt,0,PSF)

% Gamma is L*(K+M)
% F is K*T

[L, T] = size(X);

% ALS convergence choices
als_opt.MaxIterations       = 1000;
als_opt.Tolerance           = 1e-6;

if nargin>4 && ~isempty(varargin{1})
    PSF = varargin{1};
    PSF_version = true;
    M = size(PSF,1);
else
    PSF_version = false;
end

if K>0
    % Run standard IPCA 
    [Gamma_Old,s,v] = svds(double(X(:,2:T)),K);
    F_Old = NaN(K,T); F_Old(:,2:T) = s*v';
    tol         = 1;
    iter        = 0;
    while iter<=als_opt.MaxIterations && tol>als_opt.Tolerance
        [Gamma_New, F_New] = IPCA_base(Gamma_Old,X,W,Nt);
        tol       = max([ abs(Gamma_New(:)-Gamma_Old(:)) ; abs(F_New(:)-F_Old(:)) ]);
        F_Old     = F_New;
        Gamma_Old = Gamma_New;
        iter      = iter+1;
    end
    
    if PSF_version
        % Rerun IPCA with PSF
        Gamma_Old   = [Gamma_New zeros(L,M)];
        F_Old       = F_New;
        tol         = 1;
        iter        = 0;
        while iter<=als_opt.MaxIterations && tol>als_opt.Tolerance
            [Gamma_New, F_New] = IPCA_base(Gamma_Old,X,W,Nt,PSF);
            tol       = max([ abs(Gamma_New(:)-Gamma_Old(:)) ; abs(F_New(:)-F_Old(:)) ]);
            F_Old     = F_New;
            Gamma_Old = Gamma_New;
            iter      = iter+1;
        end
    end
else
    Gamma_Old   = zeros(L,M);
    F_Old       = [];
    tol         = 1;
    iter        = 0;
    while iter<=als_opt.MaxIterations && tol>als_opt.Tolerance
        [Gamma_New, F_New] = IPCA_base(Gamma_Old,X,W,Nt,PSF);
        tol       = max([ abs(Gamma_New(:)-Gamma_Old(:)) ; abs(F_New(:)-F_Old(:)) ]);
        F_Old     = F_New;
        Gamma_Old = Gamma_New;
        iter      = iter+1;
    end
end

Gamma = Gamma_New;
F = F_New;
return